%=============================================================================================
%Calibration of A and Nb to match banking data from 2014-2018
% Copy right: Yu Zhu
%=============================================================================================

clear
close all
%=====================================================================
% read all outcome
%=====================================================================

retailopt=0 % match retail to GDP 6% if value is 1 otherwise pick the best fit of the money demand curve.

k=1

foldername = '';
for i=5:99
    name = ['intermediate_results/' foldername '/output' num2str(i) '.txt'];
    if exist(name, 'file') == 2
    pars_candid(k,:) = importdata(name);
    k=k+1;
    end
    
end
foldername = 'LN87_08_final';
year_option = 0; % 1: using data from 2004-2008 0: using data from 2014-2019 2: large loan rate


figure
plot(pars_candid(:,end),pars_candid(:,1))
if retailopt==1
    [m,n] = min((pars_candid(:,6)-0.06).^2);
    pars= pars_candid(n,:);

else
    [m,n] = min(pars_candid(:,1));

pars = pars_candid(n,:);
end
%============================================================================
% define variables
%============================================================================
if year_option==1
    
%     pi_m=0.0319; % annual inflation
%     
%     i_l = 0.03826; % first percentile
%     %i_l = 0.0519061 % average loan rate
%     i_D=0.0037; % annual rate for demand deposits
%     cost=0.0265; %deposit handling cost

%% 2007-2008
     pi_m=0.03346; % annual inflation
    
    i_l = 0.03986; % first percentile
    %i_l = 0.0519061 % average loan rate
    i_D=0.004418; % annual rate for demand deposits
    cost=0.025486; %deposit handling cost
    beta = 0.96; % annual discount factor

elseif year_option ==0
    pi_m=0.01515; % annual inflation
    
    i_l = 0.0369474; % first percentile
    %i_l = 0.0519061 % average loan rate
    i_D=0.0030494; % annual rate for demand deposits
    cost=0.0202; %deposit handling cost
    beta = 0.96; % annual discount factor
else
     pi_m=0.01515; % annual inflation
    
    %i_l = 0.06; % average loan rate
    i_l = 0.0519061 % average loan rate
    i_D=0.0030494; % annual rate for demand deposits
    cost=0.0101; %deposit handling cost
    beta = 0.99; % annual discount factor

end

omega_grid = importdata('omega.txt');

alpha1 = omega_grid.data(1);
alpha2 = omega_grid.data(2);
alpha3 = omega_grid.data(3);


epsilon = 0.001; % ensures the DM utility function is 0 at 0


eta = 0.66;  %Poduction function of firms f(k)=A*k^(eta)



rq_ratio = 0.056;

reserve_rate = 0.0102; % reserve rate
pi_eff = (1+pi_m)/(1+reserve_rate)-1;
tradep = pars(end);
alpha1 = alpha1*tradep;
alpha2 = alpha2*tradep;
alpha3 = alpha3*tradep;
theta = pars(2);
sigma = pars(3);
B = pars(4);
%=================================================================
%Define functions
%=================================================================
u = @(x) ((x+epsilon).^(1-sigma)-epsilon.^(1-sigma))./(1-sigma); %DM utility
du = @(x)(x+epsilon).^(-sigma); %derivative of DM utility
d2u=@(x)(x+epsilon).^(-sigma-1)*(-sigma);
c=@(x)x;                                                         % DM cost
dc=@(x)1;%   derivative of DM cost
d2c=@(x)0;
y_star = fminsearch(@(x)(du(x)-dc(x)).^2,1);                     % y_star
g_func=@(x) (1-theta).*u(x)+theta.*c(x);                         %trading proctocol 
dg_func=@(x) (1-theta).*du(x)+theta.*dc(x);
d2g_func=@(x)(1-theta).*d2u(x)+theta.*d2c(x);
y_grid = linspace(0,y_star,1e4);
p_grid = g_func(y_grid);                                        
y_func = @(p)interp1(p_grid,y_grid,min(p,max(p_grid)));          %trading quantity 
p_func=@(x)min(x,max(p_grid));                                   %trading payment
lambda_func = @(z) max(du(y_func(z))./dg_func(y_func(z))-1,0);          % liquidity premium
D_grid = linspace(0,max(p_grid),2e3);
dlambda_func=@(z)((d2u(y_func(z)).*dg_func(y_func(z))-du(y_func(z)).*d2g_func(y_func(z)))./dg_func(y_func(z)).^3).*(z<=p_grid(end));
xi_func=@(rho)(1-rq_ratio)*max([1+rho,1/(1+pi_m)])+rq_ratio/(1+pi_m)-cost;
rho_bar = fzero(@(x)xi_func(x)-1/beta,[1e-10,2]);
rho_grid=linspace(1/(1+pi_m)-1,rho_bar,100);
%================================================
%pass variables
%================================================
otherpar.beta=beta;
otherpar.epsilon = epsilon;
otherpar.cost=cost;
otherpar.pi_eff=pi_eff;
otherpar.pi_m = pi_m;
otherpar.spread = i_l-i_D;
otherpar.i_D = i_D;
otherpar.i_l = i_l;
otherpar.alpha1=alpha1;
otherpar.alpha2=alpha2;
otherpar.alpha3= alpha3;
otherpar.rq_ratio=rq_ratio;
otherpar.eta=eta;


otherpar.lambda_func=lambda_func;
otherpar.dlambda_func=dlambda_func;

otherpar.p_bar=p_grid(end);
otherpar.rho_bar = rho_bar;





%=================================================================
% Compute inverse deman function for deposits
%=================================================================
i_m = 1/beta*(1+pi_m)-1;
psiz_grid = D_inv_demand(i_m,D_grid,otherpar);

%=======================================================================
%calibration
%=======================================================================
com_opt=1;
A = fzero(@(x)calibration_func_v3(x,otherpar,psiz_grid,com_opt),[1.0,1.6])
[drate_diff,nb,eq_final,d2c]= calibration_func_v3(A,otherpar,psiz_grid,com_opt);

pars=[sigma,B,nb,alpha1,alpha2,alpha3,beta,cost,eta,A,theta,rq_ratio,reserve_rate,year_option] 
save(['pars_post_' foldername ],'pars')
